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Abstract 

I We study in a unified fashion several quadratic vector and matrix equations with 

r—i nonnegativity hypotheses, by seeing them as special cases of the general problem Mx = 

a + b(x,x), where a and the unknown x are componentwise nonnegative vectors, M is 
a nonsingular M-matrix, and b is a bilinear map from pairs of nonnegative vectors to 
nonnegative vectors. Specific cases of this equation have been studied extensively in the 
past by several authors, and include unilateral matrix equations from queuing problems 
[Bini, Latouche, Meini, 2005], nonsymmetric algebraic Riccati equations [Guo, Laub, 2000], 
and quadratic matrix equations encountered in neutron transport theory [Lu, 2005]. 

We present a unified approach which treats the common aspects of their theoretical 
properties and basic iterative solution algorithms. This has interesting consequences: 
in some cases, we are able to derive in full generality theorems and proofs appeared in 
literature only for special cases of the problem; this broader view highlights the role of 
hypotheses such as the strict positivity of the minimal solution. In an example, we adapt an 
algorithm derived for one equation of the class to another, with computational advantage 
with respect to the existing methods. We discuss possible research lines, including the 
relationship among Newton-type methods and the cyclic reduction algorithm for unilateral 
quadratic equations. 

Keywords: quadratic vector equation, nonsymmetric algebraic Riccati equation, quasi- 
block-diagonal queue, Newton's method, functional iteration, nonnegative matrix. 
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O 1 Introduction 

O 

In this paper, we aim to study in a unified fashion several quadratic vector and matrix 
. ^ equations with nonnegativity hypotheses. Specific cases of such problems have been studied 

^ extensively in the past by several authors. For references to the single equations and results, 

we refer the reader to the following sections, in particular section 3) Many of the results 



S 



appearing here have already been proved for one or more of the single instances of the problems, 
resorting to specific characteristics of the problem. In some cases the proofs we present here 
are mere rewritings of the original proofs with a little change of notation to adapt them to 
our framework, but in some cases we extend the existing results to other problems of the 
considered class and understand the role of some key assumptions such as the positivity of the 
minimal solution. 

It is worth noting that Ortega and Rheinboldt [25, Chapter 13], in a 1970 book, treat a 
similar problem in a far more general setting, assuming only the monotonicity and operator 
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convexity of the involved operator. Since their hypotheses are far more general than those of 
our problem, the obtained results are less precise than those we are reporting here. Moreover, 
all of their proofs have to be adapted to our case, since the operator F(x) we are dealing with 
is operator concave instead of convex. 

Useful results on M-matrices In the following, A > B (resp. A > B) means Ay > B^j 
(resp. Aij > Bij) for all i, j. A real square matrix Z is said Z -matrix if < for all i ^ j. 
A Z-matrix is said an M-matrix if it can be written in the form si — P, where P > and 
s > p(P) an d p(-) denotes the spectral radius. 
We make use on the following results. 

Theorem 1. The following facts hold. 

1. If Z is a Z -matrix and there exists a vector v > such that Zv > 0, then Z is an 
M-matrix; 

2. If Z is a Z-matrix and Z > M for an M-matrix M, then Z is an M-matrix. 

3. A nonsingular Z-matrix Z is an M-matrix if and only if Z^ 1 > 0. 

4- A Z-matrix Z is a nonsingular M-matrix if and only if it has a representation as N — P, 
where N is a nonsingular M-matrix, P > and p(N~ 1 P) < 1. 

5. A Z-matrix Z is an M-matrix if and only if it has a representation as N — P, where N 
is a nonsingular M-matrix, P > and p(N~ 1 P) < 1. 

Proof. Items [l]^4] are found in [3] . We report here a self-contained proof of item |5j which was 
suggested by one of the referees of this paper. 

if Z is an M-matrix, Z = si - B with B > and p(B) < s. Then N = (s + 1)1, 
P = B + 1 is a splitting with the required properties. 

<= Let Z = N — P be a splitting with the stated properties, and let N = tl — C with 
C > 0, p{C) < t. For any e > we have piN^j^P) < 1, thus Z £ = N — j^P is a 
nonsingular M-matrix. Then Z s = tl - (C + jj^P), so t > C + j^P- Letting e ->■ + 
gives t > p(C + P). Therefore, Z = 1 1 - (C + P) is an M-matrix. □ 

Moreover, we need the following extension of item [2] 

Theorem 2. Under the hypotheses of item^ of Theorem^ if at least one of the following 
two additional conditions holds 

• M is nonsingular, 

• M is irreducible and Z ^ M , 
then Z is nonsingular. 

Proof. The results follow from the fact that the Perron value of a nonnegative matrix is an 
nondecreasing function of its entries, and a strictly increasing one if the matrix is irreducible 
®. □ 



2 



2 General problem 



We are interested in solving the equation 

Mx = a + b(x,x) (1) 

(quadratic vector equation, QVE) where M € M nxn is a nonsingular M-matrix, a, x £ M n , 
a, x > 0, and 6 is a nonnegative vector bilinear form, i.e., a map 6 : M n x M n — >■ R n such that 

1. b(v, •) and b(-,v) are linear maps for each v S M n (bilinearity) ; 

2. b(x,y) > for all x,y > (nonnegativity) . 

The map b can be represented by a tensor -B^, in the sense that b(x,y)k = Y%j=i ^ijk x iDi- 
It is easy to prove that < x < y and < z < w imply b(x, z) < b(y, w). If N is a nonsingular 
M-matrix, N~ 1 B denotes the tensor representing the map (x,y) i— > N _1 b(x,y). Note that, 
here and in the following, we do not require that b be symmetric (that is, b(x, y) = b(y, x) for 
all x, y): while in the equation only the quadratic form associated to b is used, in the solution 
algorithms there are often terms of the form b(x, y) with x ^ y. Since there are multiple ways 
to extend the quadratic form b(x,x) to a bilinear map b(x,y), this leaves more freedom in 
defining the actual solution algorithms. 

We are only interested in nonnegative solutions > 0; in the following, when referring to 
solutions of ([T| we shall always mean nonnegative solutions only. A solution x* of ([TJ is said 
minimal if x* < y for any other solution y. 

Later on, we give a necessary and sufficient condition for ([T]) to have a minimal solution. 

3 Concrete cases 

El: Markovian binary trees in [21 [TB], the equation 

x = a + b(x, x), 

with the assumption that e = (1,1,..., 1) T is a solution, arises from the study of 
Markovian binary trees. 

E2: Lu's simple equation in j2H[23], the equation 

u = u o (Pv) + e, 
v = v o (Pu) + e, 

where u, v E M. m are the unknowns, e is as above, P and P are two given nonnegative 
m x m matrices, and a o b denotes the Hadamard (component-wise) product, arises as a 
special form of a Riccati equation appearing in a neutron transport problem. By setting 
w := [u T v T ] T , the equation takes the form ([!]). 

E3: Nonsymmetric algebraic Riccati equation in |1U| . the equation 

XCX + B - AX - XD = 0, 
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where X,B £ R™^™*, C £ l m2Xmi , A £ R m i xm i, D £ R m ^ m ^ and 



M 



D -C 
-B A 



(2) 



is a nonsingular or singular irreducible M-matrix, is studied. Vectorizing everything, we 
get 

(I ® A + D T <g> I) vec(X) = vec(B) + vec(XCX), 
which is in the form Q with n = m\m2. 
E4: Unilateral quadratic matrix equation in several queuing problems [TJ, the equation 

X = A + BX + CX 2 , 

with A, B,C,X £ R mxm , A,B,C> 0, and (A + B + C)e = e, is considered. Vectorizing 
everything, we again get the same class of equations, with n = m 2 : in fact, since Be < e, 
(J — B)e > and thus / — S is an M-matrix. 

To ease the notation in the cases E3 and E4, in the following we shall set Xk = vec(Xk), 
and for E3 also m = max(mi, m-z). 

4 Minimal solution 

Existence of the minimal solution It is clear by considering the scalar case (n = 1) that 
([TJ may have no real solutions. The following additional condition allows us to prove their 
existence. We call a linear map I : R n —> R m weakly positive if l{x) > 0, l{x) ^ whenever 
x > 0, x / 0. 

Condition Al There are a weakly positive map I : R n — > R m and a vector z £ R m , z > 
such that for any x > 0, it holds that l(x) < z implies l(M~ 1 (a + b(x, x))) < z. 

We first prove a lemma on weakly positive maps, and then our existence result. 

Lemma 3. Let I : R n — > R m be weakly positive. For each z £ R m such that z > the set 
{x £ 1" : x > 0, l{x) < z} is bounded. 

Proof. For each i = 1,2, ... ,n, let ei denote the i-th vector of the canonical basis. The vector 
y = l(ei) has at least a nonzero component, let it be yj > 0. Then, for each x > such that 
l(x) < z, we must have x\ < otherwise 



which contradicts /(x) < We may repeat the argument starting from any e/~, k = 2, 3, . . . , n 
in place of e\, obtaining a corresponding bound for each of the other entries of x. □ 

Theorem 4. Equation ([TJ has at least one solution if and only if Al holds. Among its 
solutions, there is a minimal one. 
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Proof. Let us consider the iteration 

x k+ i = M^ 1 (a + b(x k ,x k )) , (3) 

starting from xq = 0. Since M is a nonsingular M-matrix, we have x\ = M~ l a > 0. It is easy 
to see by induction that x k < x k ^i~. 

Xk+i -x h = M _1 (6(x fc ,Xfc) - b(x k -i, Xk-i)) > 

since b is nonnegative. We now prove by induction that l(x k ) < 2- The base step is clear: 
1(0) = < z; the inductive step is simply Al. The sequence x k is nondecreasing and bounded 
by Lemma [3] hence it converges. Its limit X* IS £1 solution to 0. 

On the other hand, if ([I]) has a solution s, then we may choose I = I and z = s; now, 
< x < s implies M~ 1 (a + b(x, x)) < M^ 1 (a + b(s, s)) = s, thus Al is satisfied with these 
choices. 

For any solution s, we may prove by induction that x^ < s: 

s - Xk+i = M~ 1 (a + b(s, s) - a - b(xk, Xk)) > 0. 
Therefore, passing to the limit, x* < s. □ 

Taylor expansion Let F(x) := Mx — a — b(x,x). Since the equation is quadratic, the 
following expansion holds. 

F{y) = F(x) + F' x {y - x) + ^(y -x,y-x), (4) 

where F' x (w) = Mw — b(x,w) — b(w,x) is the (Frechet) derivative of F and F"(w,w) = 
—2b(w,w) < is its second (Frechet) derivative. Notice that F" is nonpositive and does not 
depend on x. 

Concrete cases We may prove Al for all the examples E1-E4. El is covered by the 
following observation. 

Lemma 5. If there is a vector y > such that F(y) > 0, then Al holds, and x* < y. 

Proof. In fact, we may take the identity map as I and y as z. Clearly < x < y implies 
M _1 (a + b(x, x)) < M~ l (a + b(y, y)) < y. It is easy to prove by induction that Xk < y. □ 

As for E2, it follows from the reasoning in |24] that a solution to the specific problem is 
u = Xq + e, v = X T q + e, where X is the solution of an equation of the form E3; therefore, 
E2 follows from E3 and Lemma [5] An explicit but rather complicate bound to the solution is 
given in |18| . 

The case E3 is treated in |10l Theorem 3.1]. Since A4 in ^ is a nonsingular or singular 
irreducible M-matrix, there are vectors v±,V2 > and u\,U2 > such that Dv\ — Cv2 = U\ 
and Av2 — Bvi = 112 ■ Let us set l(x) = Xv\ and z = V2 — A~ 1 U2- We have 

[AX k+1 + X k+l D) Vl =(X k CX k + B)v x 

<X k Cv 2 + Av 2 - u 2 < X k Dv\ + Av 2 - u 2 - 
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Since X k+ \Dvi > X k Dv\ (monotonicity of the iteration), we get X k+ \V\ < V2 — A~ 1 U2, which 
is the desired result. 

The case E4 is similar. It suffices to set l(x) = vec _1 (x)e and z = e: 

X k+1 e = (I - B)~ 1 (A + CXl)e < (J — B)~ 1 (Ae + Ce) < e, 

since (A + C)e = {I - B)e 



5 Derivative at the minimal solution 



In order to obtain a cleaner induction proof, we state the next result for a class of iterations 
slightly more general than the one that we actually need. 

Theorem 6. Let x k be the sequence generated by a fixed-point iteration of the form 

Xk+i =F{x k ) = a + Px k + b(x k , x k ), 

with a > 0, P > 0, b > 0, and suppose x k converges monotonically to x* > xq. Then p(F' x ) < 1, 
where 

K '■= P + b(x,-) +b(-,x) 
is the Frechet derivative of the iteration map. 
Proof. Let e k := — x k ; we have e k+ \ = P k e k , where 

P k :=P + b(x*,-) + b(;X k ). 

It is a classical result [20 that 



limsup y Ha?* — x k \\ < p(F^) 



(5) 



we first prove that equality holds when Po e o > 0) following the argument in |14| Theorem 3.2]. 
Since P k converges monotonically to F' , for any e > we may find an integer t such that 



p(P m ) > p«) - e, Vm> 



We have 

lim sup - x k \\ = lim sup tf\\P k -i ...Pi-.. Poe \ 

> limsup $1 jf _, fgeo 

Since Po^o > 0, PqCq > C\e for a suitable constant C\. Also, P, k 
suitable v k j > with \\v k Al = 1, so 



Pf~ l Vk,l 



for a 



limsup y |a;* — x k \\ > lim sup aC\ 



> lim sup {/Q P, fc -^ fc 



: lim sup (7 C; 



P 



fe-i 
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Since e is arbitrary, this shows that equality holds in Q. 

The case in which Poeo has some zero entries needs additional considerations. We prove 
the result by induction on the dimension n of the problem. For n = 1, either Poeo > 0, and 
thus the proof above holds, or we are in the trivial case b(u, v) = 0. Let us prove the general 
result in the case in which Poeo has some zero entries. Suppose that (up to a permutation of 
the entries) 



Poeo 



Partition conformably 



.1 :'f. 



t > 0. 



X*2 



As ei = Poeo is the error — x±, the second block row of Xk needs only one iteration to 
converge, i.e., for all k > 1, 

tk 



Xk 



X*2 



for a suitable sequence tk- Moreover, as Po > and eo > 0, for Poeo to have null components 
we need a special zero structure in P and b, namely 



P = 
Therefore, 



Pl,i Pl,2 




pL,l Pi, 2 





b(-,x* 



Ci,i Cl t 2 




F' 



Pl,l + Pl,l + Ci,l Pi, 2 + Pi, 2 + Cl )2 





/ 


•( 


tk 


5 


tk 


) 






X*2 




X*2 





and our thesis is equivalent to p{P\A + Pn,i + Ci,i) < 1- The sequence tfc, for fc > 1, converges 
monotonically to t* = x*i > ti and is generated by the fixed-point iteration 



tk+i = oi + Pi,i*fc + Pl,2^*2 + 



whose Frechet derivative at the limit point = x*i is precisely P\ t \ + B\^\ + Ci^. Therefore 
our claim holds by the inductive hypothesis. □ 

Corollary 7. By applying the theorem above to the fixed-point iteration (|3j), we obtain that 
for a QVE with x* > p(M~ 1 (b(x* : ■) + &(•, x*)) < 1 and so F' Xt is an M-matrix. 

Corollary 8. If x* > and F' Xt is irreducible or nonsingular, then by 
nonsingular M-matrix for all < x < x*, i/i,. 



Theorem 2 



P' is a 



Concrete cases For El, only the nonsingular case is of practical interest, thus the results 
are easier to prove. A strategy to reduce a problem with reducible F' x to two smaller ones is 
presented in jS] . Positivity of the solution and irreducibility are clear for E2 due to the form 
of the problem. Positivity of the solution has been proved for E3 in [TT] in the case when M 
is irreducible. Earlier versions of the results appearing in this paper |26| |2"7] contained an 
incorrect proof which failed to consider possible zero entries in Poeo- 
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6 Functional iterations 



6.1 Definition and convergence 

We may define a functional iteration for ([I]) by choosing a splitting b = b\ + 62 such that 
bi > and a splitting M = N — P such that N is an M-matrix and P > 0. We then have the 
iteration 

(N - bi(-,x k ))x k+ i = a + Pxfc + b 2 (x k , x k ). (6) 

Theorem 9. Suppose that the condition in Corollary^holds. Let xq be such that < xq < x* 
and F(xq) < (e.g. xo = 0). Then: 

1. N — bi(-,Xk) is nonsingular for all k, i.e., the iteration ^ is well-defined. 

2. Xk < Xk+\ < x*, and x& — >• x* as k ^ 00. 
5. F(xfc) < /or a// fc. 

Proof. Let J(x) := N—bi(-,x) and g(x) := a + -Px + 02(x, x). It is clear from the nonnegativity 
constraints that J is nonincreasing (i.e., x < y => J(x) > J(y)) and g is nondecreasing (i.e., 
< x < y => g(x) < g(y)). Furthermore, J(x) is a Z- matrix for all x > and J(x) > F' x . 
Under our assumptions, these results imply that J{x) is a nonsingular M-matrix for all 
< x < 7^ x* , by Corollary M 

We shall first prove by induction that x k < x* . This shows that the iteration is well-posed, 
since it implies that J{x k ) is a nonsingular M-matrix for all k. Since g(x*) = J(x*)x* < 
J(x k )x* by inductive hypothesis, ([6| implies 

J(x fc )(x* - x fc+ i) > g(x*) - o(x fc ) > 0, 

thus, since J(x k ) is ct nonsingulctr A^f-nicitrix by inductive hypothesis, — 

> 0. 

We now prove by induction that x k < x k+ \. For the base step, since we have F(xq) < 0, 
and J(xo)xo — g(xo) < 0, thus x\ = J(xo)~ 1 g(xo) > xq. For k > 1, 

J(x fc _i)(x fc+ i - x fc ) > J(xfc)x fc+ i - J(x fc _i)x fc = g(xfc) - g(xk-i) > 0. 

thus x k < Xfc+i. The sequence x k is monotonic and bounded above by x*, thus it converges. 
Let x be its limit; by passing @ to the limit, we see that x is a solution. But since x < x* 
and Xjjt IS minimal, it must be the case that x = x*. 
Finally, for each k > 1 we have 

F(xfc) = J{x k )x k - g{x k ) < J(x fc _i)x fc - g(x k -\) = 0. □ 

Theorem 10. Let f be the map defining the functional iteration ([6]), i.e. f(x k ) = J(x k )~ 1 g(x k ) = 
x k+ \. Let J, g, f be the same maps as J, g, f but for the special choice 62 = 0, P = 0. Then 
f k (x) > f k (x), i.e., the functional iteration with 62 = 0, P = has the fastest convergence 
among all those defined by Q. 

Proof. For each splitting, the functional iteration f(x) is monotonic, i.e., f{x) < f(y) whenever 
< x < y. Therefore, it suffices to prove that /(x) < /(x). We have 

x — f(x) =x — J(x) _1 g(x) = J(x)~ 1 (J(x)x — g(x)) = J(x)^ 1 F(x) 
> J(x) _1 F(x) = J(x)- 1 (J(x)x - g(x)) = x - /(x), 

which shows our claim. □ 



8 



Corollary 11. Let 

^fc+i = J{vk)' 1 gk, (7) 

where y k is a vector such that x& < y k < Xk+i, and g k a vector such that g{x k ) < 5% < g(xfc+i). 
It can be proved with the same arguments that x^+i < < x*. This implies that we can 
perform the iteration in a "Gauss-Seidel" fashion: if in some place along the computation 
an entry of x k is needed, and we have already computed the same entry of X}-+i, we can 
use that entry instead. It can be easily shown that J(x k )~ g(x k ) < J(y k )~ 1 g k , therefore the 
Gauss-Seidel version of the iteration converges faster than the original one. 

Remark 12. The iteration ^ depends on b as a bilinear form, while Equation ([I]) and its 
solution depend only on b as a quadratic form. Therefore, different choices of the bilinear form 
b lead to different functional iterations for the same equation. Since for each iterate of each 
functional iteration both x k < x* and F(x k ) < hold (thus x k is a valid starting point for 
a new functional iteration), we may safely switch between different functional iterations at 
every step. 



Concrete cases For El, the algorithm called depth in |2| is given by choosing P = 0, b\ = 0. 
The algorithm called order in the same paper is obtained in two variants with P = 0, 62 = 0, 
either on the original problem or on the one with bilinear form b(x, y) := b(y, x). The algorithm 
called thicknesses in [16] is given by performing alternately one iteration of each of the two 
above methods. 

For E2, Lu's simple iteration [24] and the algorithm NBJ in [1] can be seen as the basic 
iteration ([3]) and the iteration Q with P = 0, &2 = 0, respectively. The algorithm NBGS in 
the same paper is a Gauss-Seidel-like variant. It is shown in [15] that NBGS is twice as fast 
as NBJ in terms of asymptotic rate of convergence. 

For E3, the fixed point iterations in |14] are given by 62 = b and different choices of P. 
The iterations in [TU] are the one given by 62 = 0,P = and a Gauss-Seidel-like variant. 

For E4, the iterations in [7, chapter 6] can also be reinterpreted in our framework. 

7 Newton's method 

7.1 Definition and convergence 

We may define the Newton method for the equation ([!]) as 

K k ( x k+i ~ x k ) = -F{x k ). (8) 

Alternatively, we may write 

K k ( x k+i) = a- b(x k ,x k ). (9) 

Also notice that 

- F(x k+1 ) = b(x k+1 - x k ,x k+1 - x k ). (10) 

Theorem 13. Suppose that the condition in Corollary [§] holds. The Newton method Q 
starting from xq = is well-defined, and the generated sequence x k converges monotonically 
to x*. 
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Proof. We shall prove by induction that x k < %k+i < £*■ We have xi = M 1 a > and 
= M~ 1 a + M _1 6(x*,x*) > M~ 1 a, so the base step holds. From (10), we get 

F' Xk+i (x k+2 - Xk+i) = H x k+i - Xk,x k +i - x k ) > 0, 
thus, since F' x is a nonsingular M- matrix, x k+ 2 > %k+i- Similarly, from 

K k+1 ( x * ~ x k+2) = K x * ~ x k+ i,x* - x k+1 ), 

thus x k+ 2 < x*. The sequence x k is monotonic and bounded from above by x*, thus it 
converges; by passing ^ to the limit we see that its limit must be a solution of ([!]), hence 
x*. □ 



7.2 Concrete cases 

Newton methods for El and E2 appear respectively in [IB] and |23| . Monotonic convergence 
of the Newton method for E3 has originally been proved with the additional hypothesis 
x\ > in |14j and |10] . but this assumption was later removed in |13] . For E4, the Newton 
method is described in a more general setting in [71 122] , and can be implemented using the 
method described in [S] for the solution of the resulting Sylvester equation. However, different 
methods such as cyclic and logarithmic reduction [7] are usually preferred due to their lower 
computational cost. 



8 Modified Newton method 

Recently Hautphenne and Van Houdt [17 proposed a different version of Newton's method 
for El that has a better convergence rate than the traditional one. Their idea is to apply the 
Newton method to the equation 

G{x) = x-{M-b(-,x)y 1 a, (11) 

which is equivalent to Q. 



8.1 Theoretical properties 

Let us set for the sake of brevity R x : 
R x > F' x is nonsingular for every x < 

As for the original Newton method, it is a Z-matrix, and a nonincreasing function of x. It 
is easily seen that G' x is an M-matrix. The proof in Hautphenne and Van Houdt |T7j is of 
probabilistic nature and cannot be extended to our setting; we shall provide here a different 
one. We have 

G' x > G' Xt = R~l (M - b(-,x*) - b(R~S •)) = R x } (M - 6(.,x,) - 6(x„ •)) = 

Thus when the condition in Corollary [8] holds, G' x is a nonsingular M-matrix and thus the 
modified Newton method is well-defined. The monotonic convergence is easily proved in the 
same fashion as for the traditional method. 
The following result holds. 



= M — b(-,x). When the condition in Corollary [8] holds, 
x*. The Jacobian of G is 

= I-R x 1 b(R x 1 a,-). 
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Theorem 14 ([17]). Let Xk be the iterates of the modified Newton method and Xk those of 
the traditional Newton method, starting from x^ = = 0. Then x\~ — x\. > 0. 

The proof in Hautphenne and Van Houdt |17j can be adapted to our setting with minor 
modifications. 



8.2 Concrete cases 



Other than for El, its original setting, the modified Newton method is useful for the other 
concrete cases of quadratic vector equations. 
For E2, let us choose the bilinear map b as 



Ul 




U 2 


)- 


Ul o {PV2) 


Vl 


5 


V2 




Vl o (Pu 2 ) 



This way, it is easily seen that &(•, x) is a diagonal matrix and b(x, •) has the same structure that 
allowed a fast (with 0{n 2 ) operations per step) implementation of the traditional Newton's 
method in Bini et al. [5]. Therefore the modified Newton method can be implemented with a 
negligible overhead (0(n) ops per step on an algorithm that takes 0(n 2 ) ops per step) with 
respect to the traditional one, and increased convergence rate. 

We have performed some numerical experiments on the modified Newton method for E2; 



as can be seen in Figure 1 , the modified Newton method does indeed converge faster to the 
minimal solution, and this allows one to get better approximations to the solution with the 
same number of steps. 

For E3 and E4, the modified Newton method leads to similar equations to the traditional 
one (continuous- and discrete-time Sylvester equations), but requires additional inversions and 
products ofmxm matrices; that is, the overhead is of the same order of magnitude 0(m 3 ) of 
the cost of the Newton step. Therefore it is not clear whether the improved convergence rate 
makes up for the increase in the computational cost. 



9 Newton method and Cyclic/Logarithmic Reduction 
9.1 Recall of Logarithmic and Cyclic Reduction 

Cyclic and Logarithmic Reduction [7] are two closely related methods for solving E4, which 
have quadratic convergence and a lower computational cost than Newton's method. Both are 
based on specific properties of the problem and cannot be extended in a straightforward way 
to any quadratic vector equation. 

Logarithmic Reduction (LR) is based on the fact that if X solves 

X = B_ 1 + B 1 X 2 , 

then it can be shown with algebraic manipulations that Y = X 2 solves the equation 

Y=(I- B^Bi - BiB-i)- 1 {B_if + (/ - B^Bx - B X B„ X )- 1 (B x f Y 2 , (12) 

with the same structure. Therefore we may start from an approximation X\ = S_i to the 
solution and refine it with a term BiY, where Y is (an approximation to) the solution of (12). 
Such an approximation is computed with the same method, and refined successively by applying 
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the same method recursively. The resulting algorithm is reported here as Algorithm [TJ An 
alternative interpretation of LR [7] arises by defining the matrix- valued function / : C — >• C mxm 
as f(z) = B-i — z + B\z 2 and applying the Graeffe iteration / \-¥ f(z)f(—z), which yields a 
quadratic polynomial in z 2 with the same roots of f(z) plus some additional ones. 

Cyclic Reduction (CR) is a similar algorithm, which is connected to LR by simple algebraic 
relations (see the Bini et al. book [7] for more detail) . We shall report it here as Algorithm [2] 



9.2 Generalization attempts 

We may attempt to produce algorithms similar to LR and CR for a generic quadratic vector 
equation. Notice that we cannot look for an equation in X 2 in our vector setting, since x 2 for 
a vector x has not a clear definition — using e.g. the Hadamard (component-wise) product 
does not lead to a simple equation. Nevertheless, we may try to find an equation in b(x,x), 
which is the only quadratic expression that makes sense in our context. 

We look for an expression similar to the Graeffe iteration. If x solves = F{x) = 
Mx — a — b(x,x), then it also solves b(F(x), F(— x)) + b(F(—x), F(x)) = (notice that a 
symmetrization is needed), that is, 

b(x - M _1 a - M~H{x, x), x + M~ l a + M _1 6(x, x))+ 
b(x + M~ x a + M~ 1 b(x, x),x- M~ x a - M~ l b{x, x)) = 0. 
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Algorithm 1: Logarithmic Reduction for E4 [7] 



Input: A, B, C 

Output: minimal solution X oi X = A + BX + CX 2 

B.^il-B^A; 

B x <- (/ - B)~ l C; 

X<-B-r, 

while stopping criterion is not satisfied do 
C<- I-B-yB-x -B-1B1; 

Bx «- C- X B\\ 
X <- X + UB-x; 
U <- UBx; 
end 

return X 



Algorithm 2: Cyclic Reduction for E4 [7^ 
Input: A, B, C 

Input: minimal solution X of X = A + BX + CX 2 
R^I-B; 
S <-I-B; 
A «- A; 

while stopping criterion is not satisfied do 
5 <- R - CUT 1 A; 

x <- S'-Mo; 

Bfi-R- AR- l C - CR- l A; 
A' «- AfTM; 
C" <- CiT 1 ^ 

<- m,A',C'; 

end 

return X 
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If we set v\ = M l b(x, x) and exploit the bilinearity of &(•, •), the above equation reduces to 

- b(M ~V M^a) + (M - 6(M _1 a, •) - &(•, M^aj) v x - 6(«i, ui) = 0, (13) 

which is suitable to applying the same process again. A first approximation to x is given by 
Af _1 a; if we manage to solve (even approximately) (13), this approximation can be refined 
as x = M^a + vi. We may apply this process recursively, getting an algorithm similar to 
Logarithmic Reduction. The algorithm is reported here as Algorithm |3] It is surprising to 



Algorithm 3: A Cyclic Reduction-like formulation of Newton's method for a quadratic 
vector equation 
Input: o, Af, b 

Output: minimal solution x of Mx = a + b(x, x) 

x 0, M 4- M , a a; 

while stopping criterion is not satisfied do 

w <r- M _1 a; 

x <— x + w; 

a <— b(w, w); 

M <r- M - b(w, •) - b(-,w); 
end 

return x 



see that this algorithm turns out to be equivalent to Newton's method. In fact, it is easy to 
prove by induction the following proposition. 

Theorem 15. Let Xk be the iterates of Newton's method on ([!]) starting from xq = 0. At 
the kth iteration of the while cycle in Algorithm^ x = x^, w = x^+i — x^, M = F' Xk , 
a = -F{x k ). 



The modified Newton method discussed in section 8 can also be expressed in a form that 



looks very similar to LR/CR. We may express all the computations of step k + 1 in terms of 
Rx^b an d Rxl a only: in fact, 

Rx k+1 = Rx k ~ K; X k+ i - X k ) = R Xk [I - Rx^K; X k+1 - X k \ 

and thus 

R xl +1 Rx k = (i - Rx k b(-,x k+1 - x k )J . 

The resulting algorithm is reported here as Algorithm [4] 

The similarities between the two Newton formulations and LR are apparent. In all of 
them, only two variables {B-\ and B\, a and M, a and b) are stored and used to carry on 
the successive iteration, and some extra computations and variables are needed to extract 
from them the approximation of the solution (A, x) which is refined at each step with a new 
additive term. 

It is a natural question whether there are algebraic relations among LR and Newton 
methods, or if LR can be interpreted as an inexact Newton method (see e.g. Ortega and 
Rheinboldt [25 J, thus providing an alternative proof of its quadratic convergence. However, 
we were not able to find an explicit relation among the two classes of methods. This is mainly 
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Algorithm 4: A Cyclic Reduction-like formulation of the modified Newton method for 
a quadratic vector equation 
Input: a, M, b 

Output: minimal solution x of Mx = a + b(x, x) 

x 0, a a, b <— b, w <— 0; 

while stopping criterion is not satisfied do 

a -<— (I — &(•, w)) _1 a; 

6^(J-6_(. )U ;))-l6; 

u; <— (I — 6(a, ~ x ); 

x <— X + W] 
end 

return x 



due to the fact that the LR and CR methods are based upon the squaring X i— > X 2 , which we 
have no means to translate in our vector setting. To this regard we point out that we cannot 
invert the matrix C, since in many applications it is strongly singular. 



10 Positivity of the minimal solution 
10.1 Role of the positivity 

In many of the above theorems, the hypothesis x* > is required. Is it really necessary? 
What happens if it is not satisfied? 

In all the algorithms we have exposed, we worked with only vectors x such that < x < x* . 
Thus, if x* has some zero entry, we may safely replace the problem with a smaller one by 
projecting the problem on the subspace of all vectors that have the same zero pattern as x*: 
i.e., we may replace the problem with the one defined by 

a = ria, m = nMn T , b{x, y) = n&(n T x, n T y ), 

where II is the orthogonal projector on the subspace 



W = {x G 



for all i such that = 0}, 



(14) 



i.e. the linear operator that removes the entries known to be zero from the vectors. Performing 
the above algorithms on the reduced vectors and matrices is equivalent to performing them on 
the original versions, provided the matrices to invert are nonsingular. Notice, though, that 
both functional iterations and Newton-type algorithms may break down when the minimal 
solution is not strictly positive. For instance, consider the problem 



M = h,b 



xi 

X2 



Kxiy 2 



, X* 



For suitable choices of the parameter K, the matrices to be inverted in the functional iterations 
(excluding obviously ([3])) and Newton's methods are singular; for large values of K, none of 
them are M-matrices. However, the nonsingularity and M-matrix properties still hold for 



their restrictions to the subspace W defined in (14). It is therefore important to consider the 



positivity pattern of the minimal solution in order to get working algorithms. 
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10.2 Computing the positivity pattern 

By considering the functional iteration we may derive a method to infer the positivity 
pattern of the minimal solution in time 0(n 3 ). Let us denote by e< the t-th vector of the 
canonical basis, and es = J2sgs e s f° r an Y se t S G {1, . . . , n}. 

The main idea of the algorithm is following the iteration Xk+i = M^ 1 {a + b(xk, Xk)), 
checking at each step which entries become (strictly) positive. Since the iterates are nonde- 
creasing, once an entry becomes positive for some k it stays positive. It is possible to reduce 
substantially the number of operations needed if we follow a different strategy to perform 
these checks. We consider a set S of entries known to be positive at a certain step Xk, i.e., 
i G S if we already know that (xk)i > at some step k. At the first step of the iteration, only 
the entries i such that M a% > belong to this set. For each entry i, we check whether we 
can deduce the positiveness of more entries thanks to i and some other j G S being positive, 
using the nonzero pattern of M _1 6(-, •). As we prove formally in the following, it suffices to 
consider each i once in this process. Therefore, we consider a second set T C S of positive 
entries that have not been checked for the consequences of their positiveness, and examine 
them one after the other. 

We report the algorithm as Algorithm [5j and proceed to prove that it computes the support 
of x*. 



Algorithm 5: Compute the positivity pattern of the solution x* 
Input: a, M, b 
Output: S = {i : (z*)j > 0} 

S <— 0; // entries known to be positive 

T <— 0; // entries to check 

a' «- M -1 a; 
for i = l to n do 

if 4 > then T <- T U {i}; S <- S U {i}; 
end 

while T ^ and S / {1, 2, . . . , n} do 
t some element of T; 
T^T\{t}; 

u <(— M _1 (b(eg,et) + b(et,es))', //or only its positivity pattern 

for % G {1,... ,n} \S do 
| if m > then T <- TU {i}; S <- SU {i}; 
end 
end 

return S 



Theorem 16. The above algorithm runs in at most 0(n 3 ) operations. 

Proof. For the implementation of the sets, we shall use the simple approach to keep in memory 
two vectors S,T G {0, l} n and set to 1 the components relative to the indices in the sets. 
With this choice, insertions and membership tests are O(l), loops are easy to implement, and 
retrieving an element of the set costs at most 0{n). 

If we precompute a PLU factorization of M, each subsequent operation M~ 1 v, for v G R n , 
costs 0(n 2 ). The first for loop runs in at most O(n) operations. The body of the while loop 
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runs at most n times, since an element can be inserted into S and T no more than once (S 
never decreases). Each of its iterations costs 0(n 2 ), since evaluating b(et,es) is equivalent to 
computing the matrix-vector product between the matrix (Buj)i t j=i,..., n and eg, and similarly 
for &(e s ,ei). □ 

The fact that the algorithm computes the right set may not seem obvious at first sight. 
We prove this result by resorting to an alternative characterization of the positive entries of 
x#. For fixed a, M, b, and for a fixed i £ {1, 2, . . . , n}, we call a sequence {Sh}h=i of subsets 
of {1,2, ... ,n} positivity- showing for i if it satisfies the following properties 

i. S 1 = {h£{l,...,n}:(M- 1 a) h > 0}; 

ii. Sh C Sh+i for each h > 1; 

iii. for each t £ Sh+i \ Sh, there are r,s £ Sh such that (M B) rs t > (possibly r = s); 

iv. i £ Sn- 

Lemma 17. For each i £ {1,2, . . . , n}, (a?*)i > if and only if there exist a positivity-showing 
sequence for i. 

Proof. =^ Take i such that (x*)i > 0, and consider the iteration pi). Since x& — > x*, we 
must have (xjv)i > for sufficiently large iV. Then, we can prove that the sequence 
Sh = {i : (xfe)i > 0}, h = 1, 2, . . . , iV is positivity-showing for i. Conditions [i] and [iv] are 



clear; |n] is satisfied because the iteration is monotonic, and iii is satisfied because we 
need a nonzero summand in the right-hand side of 

(x h+l ) k = (M- l a) k + Y,{ M ~ lB )ijk{x h )i{xh)j (15) 
for the left-hand side to be positive. 

given a positivity-showing sequence, we can prove by induction on k that {xh)k > for 
each h £ Sk, where Xh are again the iterates of Q . The base step is condition [I] the 
inductive step follows from the fact that there is at least a nonzero summand in the 



right-hand side of (15) and thus the left-hand side is positive. In particular, (xat)j > 



and thus > 0. □ 

Lemma 18. The set S returned by Algorithm [5] contains i if and only if there is a positivity- 
showing sequence for i. 

Proof. =>• If at some step of the algorithm we have i £ S, then the values of S at every 
previous step of the algorithm form a positivity-showing sequence. 

<= We prove the result by induction on the length of the shortest positivity-showing 
sequence for each given i. The case iV = 1 is clear, since it must be the case that 
{M- l a)i > 0. Let us now suppose that the result is proved for all i for which the shortest 
positivity-showing sequence has length N — 1, and prove the claim for N. By condition 
there are r, s £ 5jv-i such that {M B) rs i > 0. The sequence Si, S2, ■ ■ ■ , SW-i is 
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a positivity-showing sequence of length N — 1 for all the elements of Sjv— 1> thus by 
inductive hypothesis r and s enter S (and at the same time T) at some step of the 
algorithm. If the while cycle terminates because S = {1,2, . . . ,n}, there i £ S and 
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there is nothing to prove. Otherwise, the algorithm terminates because T = 0, and thus 
r and s are removed from T at some step after being inserted. In the iteration of the 
while cycle in which either one of them is removed from T, we have m > and thus i 
enters S. □ 

The two lemmas proved above together imply that Algorithm [5] computes the correct set S. 

It is a natural question to ask whether for the cases E3 and E4 it is possible to use the 
special structure of M and b in order to develop a similar algorithm with running time 0(m 3 ), 
that is, the same as the cost per step of the basic iterations. Unfortunately, we were unable to 
go below 0(m ). It is therefore much less appealing to run this algorithm as a preprocessing 
step, since its cost is likely to outweigh the cost of the actual solution. However, we remark 
that the strict positiveness of the coefficients is usually a property of the problem rather than 
of the specific matrices involved, and can often be solved in the model phase before turning to 
the actual computations. An algorithm such as the above one would only be needed in an 
"automatic" subroutine to solve general instances of the problems E3 and E4. 



11 Other concrete cases 

In Bini et al. [6] , the matrix equation 

d 

X + Y,AiX- l Di = B-I 

i=\ 

appears, where B, Ai, Di > and the matrices B + Dj + J2i=i ^% are stochastic. The solution 
X = T — I, with T > minimal and sub-stochastic, is sought. Their paper proposes a 
functional iteration and Newton's method. By setting Y = —A -1 and multiplying both sides 
by Y, we get 

{I-B)Y = I + ^AiYDiY, 

which is again in the form 0. It is easy to see that Y is nonnegative whenever T is 
substochastic, and Y is minimal whenever T is. 

The paper considers two functional iterations and the Newton method; all these algorithms 
are expressed in terms of X instead of Y, but they essentially coincide with those exposed in 
the present paper. 



12 Research lines 

There are many open questions that could yield a better theoretical understanding of this 
class of equations or better solution algorithms. 

• Is there a way to translate to our setting the spectral theory of E4 (see e.g. Bini et al. 
chapter 3])? 

• The shift technique [7J chapter 3] is a method to transform a singular problem (i.e. one 
in which F' x is singular) of the kind E4 (or also E3, see e.g. \12\ S]) to a nonsingular 
one. Is there a way to adapt it to a generic quadratic vector equation? Is there a similar 
technique for near-to-singular problems, which are the most difficult to solve in the 
applications? 
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As we discussed in the |section 9) is there an explicit algebraic relation among Newton's 
method and Logarithmic/Cyclic Reduction, or an interpretation of the latter as an 
inexact Newton method? 

Instead of one could consider the slightly more general form 

(N - h(-,x k ) - b 3 (x k , -))x k+ i = a + Px k + b 2 (x k ,x k ), 
where 6 = 61 + 62 + 63 and M = N — P. This notation would incorporate the two 



variants of the order algorithm at the same time. We can prove as in Theorem 10 that 
P = 62 = is the best choice (among those with P, 61 , 62 , 63 > 0) , but it is not clear how 
to determine a priori the choice of 61 and 63 which gives the fastest convergence. Also, 
is there an explicit relation between the thicknesses method of Hautphenne et al. |16j 
and the symmetrized functional iteration given by 61 = 63 = ^6? 

Can this approach be generalized to the positive definite ordering on symmetric matrices 
(A > B if A — B is positive semidefinite) ? This would lead to the further unification 
of the theory of a large class of equations, including the algebraic Riccati equations 
appearing in control theory |21] . A lemma proved by Ran and Reurings [28, theorem 



2.2] could replace the first point of |Theorem 1 in an extension of the results of this paper 
to the positive definite ordering. 
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